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1 Introduction 



In recent years, much effort has been spent to design directional representation systems for 
images such as curvelets [1], ridgelets [2] and shearlets [10] and corresponding transforms 
(this list is not complete). Among these transforms, the shearlet transform stands out since 
it stems from a square-integrable group representation [4] and has the corresponding useful 
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mathematical properties. Moreover, similarly as wavelets are related to Besov spaces via atomic 
decompositions, shearlets correspond to certain function spaces, the so-called shearlet coorbit 
spaces [5]. In addition shearlets provide an optimally sparse approximation in the class of 
piecewise smooth functions with C 2 singularity curves, i.e., 

\\f - f N \\l 2 <CN- 2 (logN) 3 asiV^oo, 

where fjsr is the nonlinear shearlet approximation of a function / from this class obtained by 
taking the N largest shearlet coefficients in absolute value. 

Shearlets have been applied to a wide field of image processing tasks, see, e.g., [7, 11, 12, 17]. 
In [9] the authors showed how the directional information encoded by the shearlet transform 
can be used in image segmentation. Fig 1 illustrates the directional information in the shearlet 
coefficients. To this end, we introduced a simple discrete shearlet transform which translates 
the shearlets over the full grid at each scale and for each direction. Using the FFT this 
transform can be still realized in a fast way. This tutorial explains the details behind the 
MATLAB-implementation of the transform and shows how to apply the transform. The software 
is available for free under the GPL-license at 



http : / /www . mathematik . uni-kl . de/~haeuser/FFST/ 



In analogy with other transforms we named the software FFST - Fast Finite Shearlet Transform. 
The package provides a fast implementation of the finite (discrete) shearlet transform. 




(a) Forms with different edge ori- 
entations 



(b) Shearlet coefficients 
for a = and s = — 1 



(c) Sum 
for a 



of shearlet coefficients 



^ for all s 



Figure 1: Shearlet coefficients can detect edges with different orientations. 

This tutorial is organized as follows: In Section 2 we introduce the continuous shearlet trans- 
form and prove the properties of the involved shearlets. We follow in Section 3 the path via 
the continuous shearlet transform, its counterpart on cones and finally its discretization on the 
full grid to obtain the translation invariant discrete shearlet transform. This is different to 
other implementations as, e.g., in ShearLab 1 , see [16]. Our discrete shearlet transform can be 
efficiently computed by the fast Fourier transform (FFT). The discrete shearlets constitute a 
Parseval frame of the finite Euclidean space such that the inversion of the shearlet transform 

1 http://www. shearlab.org 
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can be simply done by applying the adjoint transform. The second part of the section covers 
the implementation and installation details and provides some performance measures. 



2 Shearlet transform 

In this section we combine some mostly well-known results from different authors. To make this 
paper self-contained and to obtain a complete documentation we also include the proofs. The 
functions where taken from [15, 14]. The construction of the shearlets is based on ideas from 
[12] and [13] . The shearlet transform and the concept of shearlets n the cone was introduced 
in [6]. 



2.1 Some functions and their properties 



To define usable shearlets we need functions with special properties. We begin with defining 
these functions and prove their necessary properties. The results will be used later. The 
following results are taken basically from [12] and [13]. 



We start by defining an auxiliary function v : 



as 



v(x) :- 



for x < 
35x 4 - 84a; 5 + 70x 6 - 20x 7 for < x < 1 

1 for x > 1. 



(1) 



This function was proposed by Y. Meyer in [15, 14]. Other choices of v are possible, in [16] 
the simpler function 



2x 2 

l-2(l-x) 2 
1 



v(x) 



for x < 
for < x < \ 
for \ < x < 1 
for x > 1 



was chosen. As we will see the useful properties of v for our purposes are its symmetry around 
(5,5) and the values at and 1 with increase in between. A plot of v is shown in Fig. 2(a). 



Next we define the function b : 



with 



'sin(f v(\oj\ - 1)) for 1 < \u\ < 2 
b{uj) := { cos(fv(||w| - 1)) for 2 < |w| < 4 



(2) 







otherwise, 



where b is symmetric, positive, real and supp6 = [—4,-1] U [1,4]. Wc further have that 
6(±2) = 1. A plot of b is shown in Fig. 2(b). 

Because of the symmetry we restrict ourselves in the following analysis to the case uj > 0. Let 
bj := b{2-i-), j € N , thus, suppfy = 2? [1,4] = [2^' +2 ] and bj(2 j+1 ) = 1. Observe that bj is 
increasing for oj 6 [2 J , 2 J+1 ] and decreasing for co G [2 J+1 , 2 J+2 ]. Obviously all these properties 
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(a) v(x) (b) solid: &(w), dashed: b(2ui) 

Figure 2: The two auxiliary functions v in (1) and b in (2) 



carry over to b 2 . These facts are illustrated in the following diagram where /* stands for the 
increasing function and \ for decreasing function. 



bj+i 



23 2i +1 
/ 1 




2 i+2 


1 \ 



2 i+3 




For j\ ^ ji the overlap between the support of b 2 ^ and b 2 2 is empty except for \ j\ —ji \ = 1. Thus, 
for b 2 and b 2 +1 we have that supp b 2 n supp b 2 +1 = [2 J+1 , 2 J+2 ]. In this interval b 2 is decreasing 

with & 2 = cos 2 (|u(^-|cj| - 1)) and b 2 +1 is increasing with 6 2 +1 = sin 2 (^v(2~^ +1 ^\u}\ - 1)). We 
get for their sum in this interval 



& 2 M + 6 2 +1 (u,) = 
Hence, we can summarize 



cos 



7T 



h 2 



l)) +s - m z (-v(2-^ 1 \oj\ 



for w < 2J+ 1 

for 2^' +1 < w < 2^ +2 

for u > 23+ 2 . 



1) 



Consequently, we have the following lemma 
Lemma 2.1. For bj defined as above, the relations 



and 



E ^) 



E ^ 



oo 

J=-l 



b 2 (2- j oj) = 1 /or |w| > 1 



sin 2 (f «(2w - 1)) /or § < |w| < 1 



/or |w| < 2 
/or | a; | > 1 



(3) 



hold true. 
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Proof. In each interval [2- ?+1 , 2- ?+2 ] only bj and j > —1, are not equal to zero. Thus, it is 
sufficient to prove that b 2 + b 2 +l = 1 in this interval. We get that 

(b 2 + b 2 +l )(u;) = b 2 {2_^) + b 2 {2_^\j) 

S 2- j [2 j + 1 , 2 j + 2 ] = [2,4] e 2- j - 1 [2 j + 1 ,2 j+2 ] = [1,2] 



= cos 2 (^v(^ ■ 2~ j oj - 1)^ + sin 2 (^(2- J - 1 w - 1)) 
(\v{2-^u - 1)) + sin 2 (\v(2~i-^ - 1)) 



= cos 2 

= 1. 

The second relation follows by straightforward computation. □ 

Recall that the Fourier transform F: L2QR 2 ) — > L2(M 2 ) and the inverse transform are defined 
by 



Ff{u>) = />):= / f(t)e~ 2 ^dt, 
.F^/M = f{t) = I f(u)e 2 ™^dw. 



Now we define the function ip\ : R — > M via its Fourier transform as 

fa(u>) := ^ 2 (2^) + 6 2 H. (4) 

Fig. 3(a) shows the function. The following theorem states an important property of ipi. 
Theorem 2.2. The above defined function ip\ has supp fa = [-4, -\] U [§, 4] and /m/^/s 

J]|^i(2- 2 ^)| 2 = 1 /or |w|>l. 

i>o 

Proof. The assumption on the support follows from the definition of b. For the sum we have 

oo 

j>0 j=0 

oo 

= ^6 2 (2- 2 ^) + 5 2 (2- 2 M, 

where -2j + 1 G {+1, -1, -3, . . .} (odd) and -2j G {0, -2, -4, . . .} (even). Thus, by Lemma 
2.1, we get 

oo 

^|^(2- 2 M| 2 = E b 2 (2-^) 
i>o j=-i 
= 1. 

□ 
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By (3) we have that 



£|&(2-* W )| 5 

j>0 



for \u\ < \ 
sin 2 (f u(2w - 1)) for § < |u| < 1 

1 for IcjI > 1. 



(5) 



Next we define a second function ip2 '■ 



again in the Fourier domain - by 



y/v(l + u) for u < 
\/v(l — uj) for u> > 0. 



(6) 



The function t/>2 is shown in Fig. 3(b). Before stating a theorem about the properties of ip2 we 
need the following two auxiliary lemmas. Recall that a function /: R — ► K is point symmetric 




(a) 



(b) V 2 



Figure 3: The functions ipi in (4) and ^2 in (6) 

with respect to (a, b) if and only if 

f(a + x)-b = -f(a-x) + b V x € M. 
With the substitution z + a — > x this is equivalent to 

f(x) + /(2a -x) = 2b Viel 
Thus, for a function symmetric to (0.5, 0.5) we have that f(x) + /(l — x) = 1. 

Lemma 2.3. The function v in (1) is symmetric with respect to (0.5,0.5), i.e., t> (x)+v (1 — x) 
1 V x G M. 



Proof. The symmetry is obvious for x < and x > 1. It remains to prove the symmetry for 
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< x < 1, we have 
v{x) + v{l — x) 

= 35x 4 - 84x 5 + 70x 6 - 20x 7 + 35(1 - x) 4 - 84(1 - xf + 70(1 - x) 6 - 20(1 - x) 7 
= 35x 4 - 84x 5 + 70x 6 - 20x 7 

+ 35 £ (f) ^ k - 84 £ + 70 £ Q - 20 £ Q 

k=0 V y k=0 V y fc=0 V y fc=0 V y 

= 1. 



Note that ip2 is axially symmetric to the y-axis. 
Lemma 2.4. T/ie function ip2 fulfills 

^1{uj - 1) + + + 1) = 1 /or |w| < 1. 



Proof. We have that 



j. 2 U(l + a;) forw<0 

2 forw>0. 



Consequently we get for < oj < 1 that 

- 1) + ^l(w) + i>l(u + 1) = u(l + oj - 1) + u(l - w) + u(l - w - 1) 

= v(oj) + v{\ — u) + u(— u;) 

=0 

= 1, 

and similarly we obtain for — 1 < to < that 

(w - 1) + ^l(w) + ^(w + 1) = u(l + a; - 1) + u(l - oj) + v(l - oj - 1) 

= v(-\oj\)+v(1 - |w|) + t;(|a;|) 
=o 

= 1. 



□ 



□ 



As can be seen in the proof the sum reduces in both cases to two (different) summands, in 
particular 



With these lemmas we can prove the following theorem 



- !) + ^iM for < a; < 1 
+ + X ) for - 1 < w < 0. 
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Theorem 2.5. The function V>2 defined in (6) fulfills 

1^2 (* + 2 M| 2 = 1 f° r M < 1, J > 0. (7) 

fc=-2J 

Proof. With S3 := 2 J o; the assertion in (7) becomes 

23 

\M k + w)| 2 = 1 for |w| < 2 j , j > 0. 

k=-2i 

For a fixed (but arbitrary) G [-2 J ',2 J '] C R we need that -1 < w* + A; < 1 for Tp 2 {oj* + k) / 
since supp ip2 = [— 1, 1]. Thus, for uj* G Z, only the summands for /c G {— — 1, — u;*, — + 1} 
do not vanish. But for k = -co* ±1 we have that uj* + k = ±1 and ^(rtl) = 0. In this case 
the entire sum reduces to one summand k = -co* such that 

J2 \Mk + ^)i 2 = \H-"* + ^)i 2 = l^2(0)| 2 = 1. 

k=-2i 

If Z and a;* > the only nonzero summands appear for k £ { , [ W *J — !}• Thus, with 
< r + := a;* — [w*J < 1, we get 

2^ 

\Mk + ^)\ 2 = |^2(-KJ +^*)| 2 + l^2(-KJ - 1 +^)| 2 = \Mr + )\ 2 + |^2(1 - r+)| 2 

k=-2i 

which is equal to 1 by Lemma 2.4. Analogously we obtain for oj* Z, < that the remaining 
nonzero summands are those for k G {[w*] , [w*] + 1}. With —1 < r~ := [a;*] + w*<0we get 

23 

£ + a/)| 2 = |V> 2 (K] + lu*)\ 2 + + 1 + ^)| 2 = |^ 2 (r-)| 2 + hfe(l + r~)| 2 . 

fc=-23 

By Lemma 2.4 and since ^(x) = ^(—x), we finally obtain 

\Mr~)\ 2 + |^2(1 + OI 2 = |^ 2 (|r-|)| 2 + |^(1 - |r"|)| 2 = 1. 

□ 

2.2 The continuous shearlet transform 

For the shearlet transform we need a scaling (or dilation) matrix A a and a s/zear matrix S s 
defined by 



8 



The shearlets ip a ,s,t emerge by dilation, shearing and translation of a function ip G L2OR ) as 
follows 

^t{x):=a-h{A- l S-\x-t)) (8) 



3 



1 



"""'HI 8 *r~T 

We assume that ip can be written as ip(ui,GJ2) = ipi((jJi)ip2C^) ■ Consequently, we obtain for 
the Fourier transform 

/ / 1 s \ \ 

3 



= aie" 2 ^^^ (awi, Va(swi + W2)) 



= ase-^^^i (awi) V2 (^a~ 5 + s 

The shearlet transform SH. 1 p(f) of a function / G L 2 (M) can now be defined as follows 

SH^(f)(a,s,t) :=(f,il> a ,s,t) 
=(/, V'o.s.t) 



f(uj)ip a:S: t(oj)doj 
=J J /(w)^i(oo;i)^2 fa"3 (~ + s )) 

The same formula can be derived by interpreting the shearlet transform as a convolution with 
the function ip aj8 (x) = tp(—A~ 1 S~ 1 x) and using the convolution theorem. 

The shearlet transform is invertible if the function ip fulfills the admissibility property 

l^(wi,w 2 )| 2 



kil 2 



duj\duj2 < 00. 



2.3 Shearlets on the cone 

Up to now we have nothing said about the support of our shearlet ip. We use band- limited 
shearlets, thus, we have compact support in the Fourier domain. In the previous section we 
assumed that ip(uji,uj2) = V'i('^'i)V ; 2(^ 2 ), where we now define ipi and ip2 as in (4) and (6) 
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respectively. With the results shown for ipi with > \ and ip2 for \uj\ < 1, i.e., \uo2\ < |wi|, 
it is natural to define the area 

& : = {(wi,W2) G R 2 : |wi| > X -, \u 2 \ < |wi|}. 

We will refer to this set as the horizontal cone (see Fig. 4). Analogously we define the vertical 
cone as 

C := {(wi.wa) G M 2 : \oj 2 \ > \, \u 2 \ > |wi|}. 
To cover all M 2 we define two more sets 

C x := {(wi,W2) G M 2 : |wi| > -, |w 2 | > = ML 

C° := {(wi,W2) G M 2 : M < l,|w 2 | < 1}, 

where C x is the "intersection" (or the seam lines) of the two cones and C° is the "low frequency" 
part. Altogether R 2 = C h U C v U C x U C° with an overlapping domain 



C U :=(-l,l) 2 \(-^) 2 - 



(9) 




Figure 4: The sets C h ,C v ,C x and C° 

Obviously the shearlet ip defined above is perfectly suited for the horizontal cone. For each set 
C K , k G {h, v, x}, we define a characteristic function XC K {^) which is equal to 1 for wGC and 
for uj G" C k . We need these characteristic functions as cut-off functions at the seam lines. We 
set 



4> h (uji,u}2) := V(wi,w 2 ) = Vi(^i)V>2 ( ^ J \e.i>- 



(10) 
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4 


UJ 2 


< ^} 




s H 


a 







For the non-dilated and non-sheared t/i h the cut-off function has no effect since the support of 
t[) h is completely contained in C h . But after the dilation and shearing we have 

SUpp^ a , Si ^|K U 2):^< |wi| < 

The question arises for which a and s this set remains a subset of the horizontal cone. For 
a > 1 we have that wi < \ is in supp^^o but not in C h . Thus, we can restrict ourselves to 
a < 1. 

Having a fixed, the second condition for supp ^> 0)S) o becomes 

— "v/ct < s + — 

— \/a — s< — <y/a — s. (11) 

W2 



Since 



< 1. 
lould 



< 1 we have for the right condition y/a — s < 1 and for the left condition — y/a — s > 
— 1, hence, we can conclude 

-l + y/a<s<l-y/a. 

For such s we have supp V> a ,s,o ^ C h , in particular the indicator function is not needed for these 
s (with respective a). One might ask for which s (depending on a) the indicator function cuts 
off only parts of the function, i.e., supp V> a ,s,o H C h / 0. We take again (11) but now we do not 
use a condition to guarantee that ^ < 1 but rather ask for a condition that allows — 

Thus, the right bound y/a — s should be larger than —1 and the left bound —y/a — s s 
be smaller than 1. Consequently, we obtain 

-l-y/a<s<l + y/a. 

Summing up, we have for |s| < 1 — y/a that supp^ a)Sj o Q C h - For 1 — y/a < \s\ < 1 + y/a parts 
of supp ip a ,s,o are also i n C V t which are cut off. For |s| > \ + y/a the whole shearlet is set to zero 
by the characteristic function. If we get back to the definition of ip a ,s,o we see that the vertical 
range is determined by ip2- By definition ^2 is axially symmetric with respect to the y-axis, in 
other words the "center" of ^2 is taken for the argument equal to zero, i.e., a~ 5 (^ + s) = 0. It 

follows that for |s| = 1 the center of tp a ,s,o is at the seam-lines. Thus, for \s\ = 1 approximately 
one half of the shearlet is cut off whereas the other part remains. For larger s more would be 
cut. Consequently, we restrict ourselves to |s| < 1. 

Analogously the shearlet for the vertical cone can be defined, where the roles of oj\ and 0J2 are 
interchanged, i.e., 

i> v (uji,u} 2 ) := 4>(uj 2 ,uji) = ^i{u 2 )^2 {/~^ * cv - ( 12 ) 

All the results from above apply to this setting. For (wi,^) G C x , i.e., \oj±\ = \ui2\, both 
definitions coincide and we define 

tp x (wi, u 2 ) := ^(wi, u 2 )xc* ■ (13) 

The shearlets iph, tpv (and V'*) are called shearlets on the cone. This concept was introduced 
in [10]. 

We have functions to cover three of the four parts of M? . The remaining part C° will be handled 
with a scaling function which is presented in the next section. 
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2.4 Scaling function 



For the center part C° (or low-frequency part) we define another set of functions. To this end, 
we need the following so-called "mother" -scaling function 



<p(u) :-- 



1 for \oj\ < \ 

cos(|f (2|w| - 1)) for \ < \oj\ < 1 
otherwise. 



The full scaling function (f> can then be defined as follows 

J f(coi) for |cji| < 1, |w2| < 



(j>(uJl,UJ2) ■■ = 



|<p(w 2 ) for \u 2 \ < 1, |^i| < 



(14) 



1 for \coi\ < 2, 1^2) < 5 

cos( |u(2|a;i| — 1)) for | < |cji| < 1, |w2| < \uj\\ 

cos(|u(2|a;2| — 1)) for \ < \lo2\ < 1, < |^2| 

otherwise. 

The decay of the scaling function (p (respectively (p) is chosen to match perfectly with the 
increase of ip\. For \u\ € [\, 1] we have by (5), that 

|^M| 2 + I^HI 2 = sin 2 (%(2M - 1)) + cos 2 (%(2|u,| - 1)) = 1. (15) 

Remark 2.6. Observe that in our setting it seems not to be useful to define the scaling function 
as a simple tensor product, namely 



$(uj) :=tp(uji)(p{uj2) 



1 

cos(fv(2|wi| - 1)) 
cos(f v{2\uj 2 \ - 1)) 



for \uj-y\ < I, \u 2 \ < \ 
for \ < \ui\ < 1, \lu2\ < 2 
for \ < \u2\ < 1, |wi| < \ 



(16) 



cos(fv(2|wi| - l))cos(fu(2|w2| - 1)) for \ < < 1, \ < \co 2 \ < 1 
otherwise. 



Fig. 5 shows the difference between both scaling functions. Obviously, the first scaling function 
aligns much better with the cones. Recently in [8] a new shearlet construction was introduced 
which is based on the scaling function in (16). We discuss the new construction in Section 3.4- 

Remark 2.7. On the other hand it is possible to rewrite the definition of the original <fi as a 
shearlet-like tensor product. We obtain a horizontal scaling function <j) h and a vertical scaling 
function <j> v as follows 

^2 \ 



<p (oji,oj2) := <f(ooi)(p - — and (/) v (uji, w 2 ) := ip(u} 2 )(p 



2lu 2 



where 



oj 2 




for |w 2 | < I 

for |c*Ji I < |w 2 | < 2|o;i| 

otherwise. 
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(a) 4>(u) 



(b) *( w ) 



Figure 5: The different scaling functions in (14) and (16). 

Thus, ¥>(^) is a continuous extension of the characteristic function of the horizontal cone 
C h . 

We set 

4>a, s ,t( x ) = M x ) = H x - *)■ 

Note that there is neither scaling nor shearing for the scaling function, only a translation. 
Thus, the index "a, s, t" from the shearlet ip reduces to "i". We further obtain 

^( w ) = e- 2 ™^^). 

The transform can be obtained similar as before, namely 

Sn <p (f)(a,s,t) = (f,<p t ). 

3 Computation of the shearlet transform 

In the following, we consider digital images as functions sampled on the grid {(^, ^r) : 
(mi, 7722) G X} with X := {(mi, 7722) : mi = 0, . . . , M — 1, m2 = 0, ... N — 1}. 

The discrete shearlet transform is basically known, but in contrast to the existing literature we 
present here a fully discrete setting. That is, we do not only discretize the involved parameters 
a, s and t but also consider only a finite number of discrete translations t. Additionally, our 
setting discretizes the translation parameter i on a rectangular grid and independent of the 
dilation and shearing parameter. See Section 3.9 for further remarks on this topic. 
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3.1 Finite discrete shear lets 



Let jo := [5 log 2 -^VJ be the number of considered scales. To obtain a discrete shearlet trans- 
form, we discretize the scaling, shear and translation parameters as 



Sj,k 

trn. 



-2j _ 



43 



3 = 0, • • . ,jo ~ 1, 



= k2~ j , -2 j <k<2 j , 



/mi m 2 \ 
V Af ' ~n) ' 



m G X. 



(17) 



With these notations our shearlets becomes ?pj t k,m(%) '■= ^aj,Sj k ,tm( x ) = ^O^a/^Sj k ( x ~ tm))- 

Observe that compared to the continuous shearlets defined in (8) we omit the factor a~ ±. In 
the Fourier domain we obtain 

where fi:={( Wl ,^): ^ = - [f J , . . . , J - 1, c 2 = - [f J , . . . , \%\ - 1} . 



\ 



(a) Shearlet in Fourier domain 



for a = I and s = — | 



(b) Same shearlet in time domain (zoomed) 



Figure 6: Shearlet in Fourier and time domain. 



By definition we have a < 1 and |s| < 1. Therefore we see that we have a cut off due to the 
cone boundaries only for \k\ = 2- ? where |s| = 1. For both cones we have for |s| = 1 two "half" 
shearlets with a gap at the seam line. None of the shearlets are defined on the seam line C x . 
To obtain "full" shearlets at the seam lines we "glue" the three parts together, thus, we define 
for I A; I =2° a sum of shearlets 
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We define the discrete shearlet transform as 
Sn(f)(n,j,k,m) := < 



(f, 4> m ) for K = 0, 
(f^j,k,m) for K e{h,v}, 
W^aJ forK=x,\k\ = 2i. 



where j = 0, . . . ,jo — 1, — 2- ? + 1 < A; < 2 J — 1 and m S I if not stated in another way. The 
shearlet transform can be efficiently realized by applying the f f t2 and its inverse if f t2 which 
compute the following discrete Fourier transforms with (9 (TV 2 log TV) arithmetic operations: 

mgl mgl 

f(m) = — - > / w e \ '^m 2 /N)/ = \ ft^e^K M + n ) m El. 

MN ^ MN ^ 

We have the Plancherel formula . 

Thus, the discrete shearlet transform can be computed for k = h as follows (observe that ip is 
real) : 

SH(f)(h,j,k,m) = </,Y^,J = J^(/,^ fcjm ) 
MN 



wen 



wen 

With ()j,k{w) := i)(4~iuji,4rikuji + 2~iu)2)f(uJ 1,(^2) this can be rewritten as 
SH(f)(h,j,k,m) = ^^hki^^'^y 

Since <7j,fc(w) G C MxN the shearlet transform can be computed as an inverse FFT of (7^, thus 

SH{f){h,j,k,m) = ifft2(&- fc ) 

= ifft2(^(4-V,4^fcwi + 2-^ 2 )/>i,w 2 )). (18) 

For the vertical cone, i.e., n = v, we obtain 

SH(f)(v,j,k,m) = ifft2(^(4-^W2,4-Jfcw2 + 2-V)/(wi,W2)) (19) 
and for the seam line part with |fc| = 2 J we use the "glued" shear lets and obtain 

SU{f)^ v {j,k,m) = ifft2(^ xw (4-Jwi,4-^a;i + 2-^W2)/(wi,W2)). (20) 
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Finally for the low-pass with 50(^1^2) := <f>(wi, ^>2)f(^i, ^2) the transform can be obtained 
similar as before, namely 

sn 4> {f){m) = {f, 4> m ) 



— }_^e 'W"^( Wl , W2 )/( W i,w 2 ) 
l[N E e +2m(W ' l ^ j ^(^i,o;2)/>i,a;2) 



J_^+™<<£#)>, 



MN 



1 +27ri(^,( m i / *n>-~ / 



= if ft2(g ) 

= ifft2(^(wi,W2)/(wi,W2)). (21) 



The complete shearlet transform is the combination of (18) to (21). We summarize 



ifft2(<f)(oJi,oj2)f(ui, CO2)) for k = 

ifft2(^(4-J'wi,4- J 'A!Wi+2-Jw 2 )/(wi,t«; 2 )) for « = h, \k\ < V - 1 
o H{ i)(k, j, k,m) = < ~, „ ~ , , . (22) 

' iff 1 2 (^(4-^2, 4-Jfcw 2 + 2-^i)/(wi,cj 2 )) for K = v, \k\ < V - 1 V ; 

,ifft2(^ /lx, '(4-Jcji,4^fccji +2-^ 2 )/(wi,w 2 )) for k 7^ 0, |fc| = 2^'. 



3.2 A discrete shearlet frame 

In view of the inverse shearlet transform we prove that our discrete shearlets constitute a 
Parseval frame of the finite Euclidean space L2 (X) . Recall that for a Hilbert space T~L a sequence 
{uj : j G J} is a frame if and only if constants < A < B < 00 exists such that 

A\\f\\ 2 <^2\(f^ 3 )\ 2 <B\\f\\ 2 for all f EH. 

The frame is called tight if A = B and a Parseval frame if A = B = 1. Thus, for Parseval 
frames we have that 

II/I| 2 = EK/'^>| 2 for all /G^ 
which is equivalent to the reconstruction formula 

f = ^2(f,u j )u j for all / G U. 

Further details on frames can be found in [3] and [14]. In the n-dimensional Euclidean space 
we can arrange the frame elements Uj, j = 1, . . . , n > n as rows of a matrix U. Then we have 
indeed a frame if U has full rank and a Parseval frame if and only if U T U = I n . Note that 



16 



UU T = I n is only true if the frame is an orthonormal basis. The Parseval frame transform and 
its inverse read 

{{f,u j )% 1 = Uf and f = U T ((f, Uj })] =1 . 
By the following theorem our shearlets provide such a convenient system. 

Theorem 3.1. The discrete shearlet system 

{4 m H: j =0,...,j -l,-2 J + l<K2 J -l,meI} 
U {iplk,m(u) ■ 3 = °> • • • . Jo - h -2 j + 1 < k < V - 1, m G 1} 
U : i = 0, . . . , jo - 1, |fc| =2\mel} 

U {</> m (w) :m£l} 

provides a Parseval frame for L 2 (T). 

Proof. We have to show that 

j'o-i io-l 

= E E E Eia^)i 2 +E E Ei(/-tti 2 +Ei(/.wi J = :C - 

fce{M} 3=0 fc=- 2 i+lmeX j=0 k=±2i mel mel 



Since \\f\\ 2 F = jy^vll/lli? (Parsevals formula, ||-|| denotes the Frobenius norm, i.e., 
||vec(-) H2) it is sufficier 

By (18) we know that 



I vec ( • ) 1 1 2 ) it is sufficient to show that C is equal to j^ll/llf?- 



We further obtain 



wen 



E \(f^km)\ 2 = E \9j,k(m)\ 2 = \\9j,k\\ 2 F . 

mgl m£l 



Consequently, with Parsevals formula 



toll! - j^WkkWl - E l&.*( w )l 2 

wen 

= WnT. I^(4-V,4-^o;i + 2-^ 2 )/> 1 ,o; 2 )| 2 

= /jjjy E 1^(4-^1,4-^1 + 2-^ 2 )| 2 |/>i^ 2 )| 2 . 
wen 

Analogously we obtain for the vertical part 

E K/.^fc,m>| 2 = ^ E l^(4-^2,4-^ W2 + 2-V)| 2 |/K^2)| 2 . 
mgl wen 
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Using these results we can conclude for the seam-line part 

E K/>^Ol 2 = Wn ^ I^W"''*"! + 2-^ 2 )| 2 |/>i^ 2 )| 2 x ch 

+ ^ E W^2,4- J fc W2 + 2- i Wl )| 2 |/>i,W2)|V 
wen 

^- ^ 1^(4-^,4-^^ + 2-^ 2 )\ 2 \f(u 1 ,u J2 )\ 2 x c ,. 



+ MiV 

wen 



For the remaining low-pass part we get similarly 

^E IE £^)/(«)l s 



MiV 



MiV 

m£l wen 



1 Y^V^ 27ri (^.( mi/ /«))l/ u"/ Mi 

^EIE e Vm 2 /Ar; V(^l,^2)/(wi,W 2 )| 



MiV 

men wen 



with <?oM := <K^i, u} 2 )f(ui, u 2 ) 

IXfAfE 6 U2/Ar %H| 



'MiV 
mei wen 



|2 _ ||„ 1 1 2 1 



= E \90(m)\~ = \\90\\f = Mff\\9o\ 

= E I<AK^2)/(W1,W 2 )| 2 

wen 

= ^Ei^'^)i 2 i/>i^ 2 )i 2 . 

wen 
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Lets put the pieces together: 

j'o-i 23-1 j -l 

c = E E E Ei</.«i a + EEEi(/.0 2 + Ei(/.wi 2 

«e{M} j=0 fc=-23+lmeX j=0 fc=±2J ' mel mel 

JO-1 23-1 

= E E ^Ei^( 4 "^' 4_J '^ + 2 "^)i 2 i/>i^2)i 2 

3=0 k=-2i+l wen 
jo-1 2^-1 

|2 



JU ^ J- 

+ E E ^Ei^ 4 ~^4-^ 2 + 2-^)1 2 i/>i^ 2 ) 

3=0 fc=-2i+l wen 

+ E E (^Ei^ 4 "^i' 4 ^+ 2 ^)iV>i^2)iv 

j=o fc=±2J ^ wen 

+ ^ E 1^(4-^,4-^^2 + 2-V)| 2 |/>i, W2 )|V 
wen 

+ jj^ E l^( 4 " J wi, 4- J 'fcwi + 2-^ 2 )| 2 |/>i^2)| 2 Xcx j 
wen / 

+ ^El^'^)| 2 |/>i,o; 2 )| 2 . 



wen 

We can group the sums by the different sets and obtain 
jo-i 2J 

1AA ^ 



-1 JKJ 

c = ^]vEE E hfa-w-'k* + 2-^)i 2 i/( W i >W2 )iv 

wen j=o fc=-23 

1 JO-l 23 

+ ^EE E 1^(4-^,4-^ + 2-^)| 2 |/(^,W2)| 2 xc- 
wen j=o fe=-23 

1 jW 1 

+ ^EE |/>i^2)| 2 |^(4-V,0)| 2 x C x + — ^|^i, W2 )| 2 |/>i,a; 2 )| 



MN t—' ' yi V v v - MiV 

wen j=o _^ wen 

Using the definition of ?/> in (10) (or (12) and (13), respectively), we can conclude 

C =MN E l/>i^ 2 )| 2 ^|^i(4-V)| 2 £ \M^ + k)\ 2 
weC 1 j=o fc=-23 1 

/s ' 



=1 for |wi|>l =1 (see Theorem 2.5) 

(see Theorem 2.2) 

1 - V 

+ MN E l/>i,^)| 2 E l^( 4_J ^)| 2 E l^(2^ + fc)| 2 
wec« j=o fc=-23 2 



=1 for |w 2 |>l =1 

+ ^E l/V,^)| 2 + ^E i^wi.wa)! 2 !/^,^)! 5 

r,jf=D. 

1 for w e 



wec x wen r . 1 2 
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With the properties of ipi and ip2 (see Theorems 2.2 and 2.5) we obtain two sums, one for the 
overlapping domain C u (see (9)) and one for the remaining part 

-. _ /j'o-i io-i 

C =MN E I/>1^2)| 2 +E |/(W1,W2)| 2 El^( 4 "^)| 2+ E l^l(4-^2)| 2 + \k^2) 

ujen\c D cuec a \j=o j=o 



where we can split up the second sum as 

c = wn E l/(^^)l 2 



MiV 



+ ^ E l/>i,^)| 2 sin 2 (^(2^1-1))+^ l/>i^2)| 2 sin 2 (^(2|o; 2 |-l)) 

wec h nc n " uiec v nc a 

+ WN E l/(^i^2)| 2 cos 2 (%(2|a; 1 |-l)) + ^ £ |/>i,u> 2 )| 2 cos 2 (^(2|u, 2 | - 1)) 



using the overlap (see (15)) we can continue 



MiV 



+ ^ E l/>i^2)| 2 (sin 2 (%(2| W1 | - 1)) +cos 2 (| V (2|wi| - F 



=1 (see (15)) 



+ ^V E l/>i^2)| 2 (sin 2 (%(2|a; 2 |-l))+cos 2 (^(2|a; 2 |-l))). 
wec^nc s v ' 

=1 (see (15)) 

Finally, we obtain 

c = m77 E l/( w i' w 2)| 2 = j^f\\f\\ 2 F = II/IIf- 

wen 



□ 



3.3 Inversion of the shearlet transform 

Having the discrete Parseval frame the inversion of the shearlet transform is straightforward: 
multiply each coefficient with the respective shearlet and sum over all involved parameters. As 
an inversion formula we obtain 

JO-l 2^-1 jo-l 

f= E E E E</'V^j^ fe) ™ + E E E(/^S:j^ + E(/^-)^- 

ne{h,v} j=0 k=-2i+lm£l j=0 k=±2i m£l m£l 

The actual computation of / from given coefficients c(k, j, k, m) is done in the Fourier domain. 
Due to the linearity of the Fourier transform we get 

j"o-l 2^-1 jo-l _ _ 

/= E E E E</^,jfe,™ + E E E</> + E</>^- 
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We take a closer look at the part for the horizontal cone where we have 

jo-i 2J-1 

j=0 k=-23+l mel 

= ^ c(/i,j,fc,m)e na,,( -2/Jv^^(4-^i,4^wi + 2-^ 2 ). 

j=0 k=-2i+l mel 

The inner sum can be interpreted as a two-dimensional discrete Fourier transform and can be 
computed with a FFT and we can write 

jo-i 2 j -1 

/Mxc = 2 1^ fft2(c(/i,j,/c,-))(wi,w 2 )V'(4-V,4 i A;a;i + 2^w 2 ). 

j=0 fc=-2J+l 

Hence, / can be computed by simple multiplications of the Fourier-transformed shearlet co- 
efficients with the dilated and sheared spectra of tp and afterwards summing over all "parts" , 
scales j and all shears k, respectively. In detail we have 

/(W1,W2) =fft2(c(0,-))^(wi,W2) 

jo-1 2J -1 

+ 2 iit2(c(h,j,k,-))^( y A- j u 1 ,4~ j kuj 1 + 2' j uj2) 

j=0 k=-2i+l 

jo-i 2^-i ^3) 

j=0 fc=_23+l 
jo-1 

+ ^ fft2(c(/i x w,j,fc,-))^(4 _:? wi,4~ J 'A;a;i + 2~ j o;2)- 

Finally we get / itself by an iFFT of / 

/ = ifft2(/). 

3.4 Smooth shearlets 

In many theoretical and some practical purposes one needs smooth shearlets in the Fourier 
domain because such shearlets provide well-localized shearlets in time domain. Recently, in 
[8] a new shearlet construction was proposed that provides smooth shearlets for all scales a 
and respective shears s. Our shearlets are smooth for all scales and for all shears |s| / 1. 
The "diagonal" shearlets ip hxv are continuous by construction but they are not smooth. Fig 
7(a) illustrates this. Obviously our construction is not smooth in points on the diagonal. The 
new construction circumvents this with "round" corners. To this end, we get back to the 
two different scaling functions which we discussed in Section 2.4. While we chose the scaling 
function matching our cone-construction the new construction is based on the tensor-product 
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(a) diagonal shearlet in our construction 



(b) diagonal shearlet in the new construction 



Figure 7: Diagonal shearlets in our construction and in the new, smooth construction (Fourier 
domain) 



scaling function <E>(u;) = (p(oji)<p(u)2) ■ We present the basic steps in the construction of [8] 
transferred to our setting. In fact, we only need to modify the function ipi. We set 



$iH := y $ 2 (2~ 2 wi, 2- 2 u 2 ) - $ 2 (wi, w 2 ). (24) 
Clearly, ^(w) fulfills £\> *i(2" 2j "w) = 1 for all u € Q \ [-1, l] 2 . We further have that 



$ 2 (w) + ^ ^>l(2' 2j oj) = 1 for w G O, 
i>o 

i.e., this setting provides also a Parseval frame. Fig. 8 shows Note that is supported 
in the Cartesian corona [— 4, 4] 2 \ [—\, \} 2 - The full shearlet \P is similar as before: 



^(W1,W 2 ) = *l(wi,W 2 )^2 



U>2 
Ml 



(25) 



The construction of the horizontal, vertical and "diagonal" shearlets is the same as before, 
besides that the diagonal shearlets are smooth now, see Fig. 7(b). 

Before we examine the smoothness of the diagonal shearlets we discuss the differentiability of 
the remaining shearlets. Due to the construction we only need to analyze the functions ipi and 
ip2- We have 

for | Wl | < \ 
sin(|v(2|a;i| - 1)) for \ < < 1 

1 for 1 < |wi| < 2 
cos(ff(||(Ji| - 1)) for2<|^i|<4 
for \uji\ > 4 



^i(wi) = ^ 2 (2wi) + 6 2 K) 
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Figure 8: The new function in (24) 



and with straight forward differentiation 




for \uji\ < j 

7rv(2|wi| - l)u'(2|wi| - 1) cos(f v (2|wi| - 1)) for § < < 1 
for 1 < |cji| < 2 

-f - iy(^|wi| - 1) sin(f u(||wi| - 1)) for 2 < | Wl | < 4 

for Iwil > 4. 



The derivative is continuous if and only if the values at the critical points |, 1, 2, 4 coincide (for 
symmetry reasons we can restrict ourselves to the positive range). We have that v(2 ■ \ — 1) = 
v(Q) = (even v'{0) = 0) and v'(2 ■ 1 - 1) = v'{l) = and further v{\ ■ 2 - 1) = u(0) = 
and ■ 4 — 1) = t>'(l) = 0. Consequently, ^ is continuous and in particular G C 1 . By 
induction we see that tjj^ G C n if and only if v^ n \0) = and v^ n \l) = 0, n > 1. 

For our u in (1) we have z/ 3 )(l) = but t)' 4 '(l) / 0, i.e., ipi G C 3 . 

Similarly, we obtain for ^2 that 



9^2 

doji 



(^1,^2) 



^' 2^(1+^) 



for ^ < 
for ^ > 0, 



where we see that we need v'(0) = for the derivative to exist. Thus, the shearlet V> is C n if 
v ( n )(0) = i;( n )(l) = 0. This is also valid for the dilated and sheared shearlet 4>j k m (and $j km ) 
for I A; I / 2 J . We take a closer look at the diagonal shearlet for k = — 2 J where we have 





-2i,m 




for 




G C h 




-2J,m 




for 


a; 


G C v 




-2J,m 




for 


Co' 


GC X 
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Naturally, m is smooth for u G C h and oj G C. Additionally, m (w) is continuous at 



the seam lines, but not differentiable there since we have for the partial derivatives of m 
and ih v oi (lo) that 



( W ) = yi 2 (2^ - 1))e -¥(™+™) . 2 - 2 ^(2- 2 V) 



+ ^ 1 (2- 2 ^ 1 )e-^(™ + ™)(-2^)^(2^(^-l)) 
+ ^i(2- 2 ^i)^ 2 (2- ? '(— - l))(- 2 ^ mi ) e -¥(^ m i+™) 

TV 



and 



d<4> v 



'i^(u J )=M2' 2 ^ 2 ) e -^™+™)(-)^(2^ - 1)) 

+ 4> 2 (2 j (— - i))(-- mi ) e -f^ ra i+«) 

UJ 2 N 

For wi = w 2 we obtain 



(wi.wi) = e -^(»i+^) • (^(0)2- 2 ^(2- 2 ^ 1 ) 

=i 

- ^(2- 2 V)(^) ^(0) -^ 1 (2- 2 V)^2(0)(^m 1 )) 



=i 



• (2-^^(2-^) - (^m 1 )V'i(2- 2 ^ 1 )) 



<9u;i 

and 



^=^(^1,0,1) = e-¥-(-+-) . (^ 1 (2- 2 ^ 1 )(-) ^(0) -^i(2- 2 ^i)V' 2 (0)(^m 1 )) 



,2m 



= e -^ Wl (m 1+ m 2 ) . ^_ ( ^ mi) ^ l(2 -?, Wl) y 



Obviously, both derivatives do not coincide, thus, our shearlet construction is not smooth for 
the diagonal shearlets. If we get back to the new construction we obtain for the both partial 
derivatives 

r 2J ' m H = 2- 2 ^(2- 2 M^ 2 (2 J (^ - l)) e -^(™ + ™) 

OUJ\ OUJ\ COi 

- 2^* 1 (2- 2 M^ 2 -(2^(^ 2 - - i) )e -¥K-i+™) 

- ^mi^ 1 (2- 2 ^)V» 2 (2^(^ 2 - i)) e -¥(^ m i+™) 
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and 



r 2J ' m H = 2- 2 ^(2-^u)M2 J (- ~ l)) e -¥(™+™) 

OUJi OUJ\ UJ 2 



V W2 9wi 1^2 

2 ^m 1 i> 1 (2- 2 ^)M^(— ~ l)) e -¥(™+- 2 ™ 2 ). 

iV £J2 



With oji = u)2 we obtain further 



and 



=1 

-2^^ 1 (2- 2 ^ 1 ,2- 2 ^ 1 )^(0)e-^( mi+m2 ) 

=o 

- ^™mi*i(2- 2 ^i, 2~ 2 V) ^2(0) e -¥^( m i+ m 2) 

=i 



j^V .afr) = 2- 2 ^(2^ 1 ,2- 2 ^ 1 )^ 2 (0)e-^(™ 1 +™ 2 ) 

=i 



+ ( ^)^ l( 2-2^ 1)2 - 2 V)^(0) e -^ l(mi+m2) 

=0 

- ^mi*i(2- 2i wi, 2~ 2 V) ^j(O) e -¥^( m i+ m 2). 

=i 

It can be easily seen that both derivatives coincide if and only if g^(0) = since then the 
second term vanishes. The same result is obtained for partial derivative with respect to u)2- 
Consequently, the new construction is smooth everywhere. 

Remark 3.2. As we have seen the smoothness of the shearlets depends strongly on the smooth- 
ness of the function v. The v we have used was constructed to provide shearlets in C 3 . The 
first three derivatives at and 1 should be equal to zero, i.e., v'(x) = c • x 3 (x — l) 3 . With 
v(l) = 1 and straight forward integration one obtain c = —140 and the function v as in (1). 

Higher grades of smoothness can easily be obtained by creating a new v by setting v'(x) = 
c ■ x k (x — l) k . These shearlets would be in C k . 

Note that due to our discretization t = m we have a unique handling of both the horizontal 
and the vertical cone and do not have to make any adjustments for the diagonal shearlets, in 
contrast to the discretization t = A aj S Sjk m where you have different discretizations for t in 
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the horizontal and in the vertical cone. Consequently, some adjustments have to be made for 
the diagonal shearlets. 

Smooth shearlets are well- located in time. To show the difference we present in Fig. 9(a) the 
"old" shear let in the time domain and in comparison in Fig. 9(b) the new construction in the 
time domain. The non-smooth construction is slightly worse located. The shearlet coefficients 



Figure 9: diagonal shearlets in our construction and in the new, smooth construction (Time 
domain) 

of, e.g., a diagonal line only show marginal differences such that for practical applications it is 
irrelevant which construction is used. 

3.5 Implementation details 

The implementation of the shearlet transform follows very closely the details described here. 
As we see in (22) and (23) for both directions of the transform the spectra of tp and <j> are 
needed for all scales j and all shears s on "all" sets. We precompute these spectra to use them 
for both directions. 

Up to now our implementation only supports quadratic images, i.e., M = N (see Remark 3.4 
for a short discussion on this topic). 

3.5.1 Indexing 

To reduce the number of parameters we introduce one index i which replaces the parameters k, 
j and k. We set i = 1 for the low-pass part. We continue with the lowest frequency band, i.e., 
j = 0. The different "cones" and shear parameters represent the different "directions" of the 
shearlet. For illustration we reduce the shearlet to a line which is rotated counter-clockwise 
around the center and assign the index i accordingly. In each frequency band we start in the 




(a) diagonal shearlet in our construction, time 
domain 



(b) diagonal shearlet in the new construction, 
time domain 
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horizontal position, i.e., k = h and k = 0, and increase i by one. For each k = —1, . . . , —2 3 + 1 
we continue increasing the index by one. The line is now almost in a 45° angle (or a line with 
slope 1). The next index is assigned to the combined shear let u h x v" at the seam line which 
covers the "diagonal" for k = —2 3 . We continue in the vertical cone for k = — 2 J ' + 1, . . . , 2 J ' — 1. 
Next is again the combined shearlet for k = 2 3 . With decreasing shear, i.e., k = 2 3 \ . . . , 1, we 
finish the indexing for this frequency-band and continue with the next one. 

Summarizing we have always one index for the low-pass part. In each frequency band we have 2 
indices (or shear lets) for the diagonals (k = ±2 3 ). In each cone we have 1 + 2 • (2* - 1) = 2 3+1 - 1 
shearlets. For the scale j we have 2 • (2 3+l — 1) + 2 = 2 3+2 shearlets. The following table lists 
the number of shearlets for each j. 



low-pass 


3 = 


J = l 


i = 2 




1 


4 


8 


16 





With a maximum scale jo — 1 the number of all indices rj is 

jo-i io-i 
7] = 1 + 2j+2 = 1 + 4 ^ 2 J ' = 1 + 4 • (2 30 - 1) = 2 jo+2 - 3. (26) 

j=0 j=0 

For each index the spectrum is computed on a grid of size N x N. We store all indices in a 
three-dimensional matrix of size N x N x r/. The first both components refer to the o>2 and oj\ 
coordinates and the third component is the respective index. Consequently, an image / of size 
N x N is oversampled to an image of size N x N x r/. In particular we have an oversampling 
factor of rj. The following table lists r\ for jo = 1, . . . , 4. 



Jo 


1 


2 


3 


4 




V 


5 


13 


29 


61 





Note that jo is the number of scales, the highest scale parameter j is always jo — 1, i.e., we 
have the scale parameters 0, . . . , jo — 1. The function helper/shearletScaleShear provides 
various possibilities to compute the index i from a and s or from j and k and vice versa. See 
documentation inside the file for more information. 

3.5.2 Computation of spectra 

We compute the spectra ipj,k,m as discrete versions of the continuous functions, i.e., we compute 
the values on a finite discrete lattice S of size N x N. With the functions defined as in (4) and 
(6) we may not take S = 17 as this would destroy the frame property. The question is how to 
choose X (such that E G [— X, X) NxN ) or the distance A between to grid points, respectively. 

The use of Theorem 2.2 in the proof of Theorem 3.1 was not completely correct. We have 
that supp^iM = [\A] = [2~\2 2 ] and tpi = 1 for uj £ [1,2] = [2°,2 1 ]. For the scaled 
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version we further have supp^i(2- 2 M = [2 2 - 7-1 , 2 2 -?'+ 2 ] and = 1 for uj £ [2^,2^' +1 ]. 

Consequently, we can conclude 



io-i 
j"=o 







for IcjI < h 



for \ < \u\ < 1 



sin 2 (f v(2w - 1)) 
1 for 1 < |w| < 2 2 ^'°- 1 ) +1 

2 (| v (2-2jo-i w _ i)) f or 2 2(?o-i)+i < | w | < 2 20o-i)+2 

for M > 2 2 ^'-i)+2. 



cos 




We see that the sum is equal to 1 for a wide range of a;. As we have seen the part for \u\ < 1 
where the sum increases from to 1 matches with the decreasing part of the scaling function 
(see (15)). But we also have a decay for \oj\ > 2 2 ( J ' 0-1 ) +1 where there is no compensation to 1 
since there are no higher scales. Keeping this decay would violate the frame property. 

Consequently, X must be less or equal than 2 2 ( J0_1 ) +1 = 2 2 -? 0-1 which implies the decay to be 
"outside" the image. We set X = 2 2j0_1 . Additionally, if we get back to the lowest scale it is 
reasonable to set A < 1 since otherwise there would be to less grid points in supp^i- 

To compute the grid and the spectra we assume that A = 2n + 1 is odd. We then have a 
symmetric grid around 0, hence, we have n grid points in the negative range and n grid points 
in the positive range and one grid point at 0. If N is given even we increase it by 1. After 
computing grid and spectra we neglect the last row and column to retain the original image 
size. Having n = grid points for the positive range and the maximal distance between 
two grid points A = 1 we get 



X = 2 2jo ~ 1 = 



N ■ 



j = -log 2 (A-l). 



Since jo is the (scalar!) number of scales, we set for the number of scales (as used above) 
jo := [^log 2 {N)\. In the following table we list the number of scales for all image sizes 
N = 4..., 1024 



N 


4,..., 15 


16,..., 63 


64, ... , 255 


256, . . . , 1023 


1024 


jo 


1 


2 


3 


4 


5 



With jo fixed we can compute the second parameter A. As we have seen the highest value in 
the grid should be X = 2 2j ' 0-1 . For an odd the grid ranges from [— X, X] and for an even 
grid we have the range [—X, X) = [—X, X — A]. We assume again an odd N, thus, the interval 
[—X, X] should be divided in N grid points including the bounds —X and X what leads to 
N — 1 subintervals and 

2 • X 2 • 2 2 -?"- 1 2 2 -? 

A 



N - 1 



Where A = 1 if N = 2 2 ^ 



1 and A > |, i.e.. 



N 
i 



1 



N - 1 



4 < A < 1. Thus, for the same number of scales 
we obtain a better resolution with increasing image size. 



It seems a little awkward to discretize / and ip on different lattices. However, with this 
auxiliary construction the definition and properties of the shear let ip is much more convenient. 
Additionally the shearlets are now independent of the parameter A (or other grid properties). 
Anyway, to circumvent the imperfection with two lattices we can now formerly also discretize 
V>(Au>) on Q instead of *p(t>j) on E and obtain the same spectra. 
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Remark 3.3. The spectra depend strongly on the image size. In particular, if we reduce the 
image size a bit but still have the same number of scales the resolution of the different frequency 
bands varies. This may be not wanted for comparison issues. To circumvent this one could 
chose the grid according to the highest image size (for the respective number of scales) and drop 
the boundaries for smaller images. This would however lead to a very small high frequency band 
for smaller images. In our current implementation this frequency band is as large as possible. 

Remark 3.4. The theoretical results shown are valid for both square images and rectangular 
images, i.e., we have I and£l of size MxN. However, our implementation only supports square 
images. For the extension to rectangular images some questions remain. The first question 
is what the "diagonal" in a discrete rectangular setting is, i.e., the "diagonal" shearlets have 
to be handled carefully. More tricky is the question how to handle the different sizes in the 
both directions with respect to the size of the grid especially if M <C N (or M S> N ) and 
the number of scales. Possible are an equispaced grid in both directions and thus less scales 
in one direction or the same number of scales but a non- equispaced grid what would lead to 
rectangular frequency bands. Depending on the application both seems useful. We hope to 
implement rectangular shearlets in a future version of the software. 

3.6 Short documentation 

Every file contained in the package is commented, see there for details on the arguments and 
return values. Thus, we only want to comment on the two important functions. 

The transform for an image AG £, NxN is called with the following command 

[ST,Psi] = shearletTransf ormSpect (A , numOf Scales , shearlet_spect 
, shearlet_arg) 



where numOf scales, shearlet_spect and shearlet_arg are optional arguments. If only A is 
given the number of scales jo is computed from the size of A, i.e., jo = [3 log 2 (./V)J and the 
shearlet with (4) and (6) is used. On the other hand numOf Scales can be used two-fold. If 
given as a scalar value it simply states the number of scales to consider. On the other hand 
we can provide precomputed shearlet spectra which are then used for the computation of the 
transform. Observe that the shearlet spectra only depend on the size of the image and the 
number of scales, thus, they are cached and reused if the function is called with an image 
of same size again. The variable ST contains the shearlet coefficients as a three-dimensional 
matrix of size N x N x r] with the third dimension ordered as described in section 3.5.1. Psi 
is of same size and contains the respective shearlet spectra ipjko- 

With the parameters shearlet_spect and shearlet_arg other shearlets can be used to com- 
pute the spectra. Included in the software is the 'meyerShearletSpect' as default shearlet 
(based on (4) and (6)) and 'meyerNewShearletSpect' for the new smooth construction (see 
(25)). The parameter shearlet_arg is not used in both cases. 

The application of other shearlet spectra is very straight forward. One the one hand one 
can simply compute them externally in the matrix Psi and provide them as the parameter 
numOf Scales. On the other hand it is possible to provide an own function 'myShearletSpect' 
(with arbitrary name) with the function head 
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Psi = myShearletSpect (x , y , a , s , shearlet_arg) 



that computes the spectrum Psi for given (meshgrids) x and y for scalar scale a and shear 
s and (optional) parameter shearlet_arg. For shearlet_arg=' scaling' it should return 
the scaling function. To obtain a reasonable transform the shearlet should provide a Parseval 
frame. To check this just compute (and plot) sum(abs(Psi) . ~2,3)-l. The values should be 
close to zero (see Fig. 11(a)). Call the shearlet transform with the new shearlet spectrum by 
setting the variable sherlet_spect to myShearletSpect or whatever you chose as the name 
of your shearlet function. 

The inverse transform is called with the command 

A = inverseShearletTransf ormSpect (ST , Psi , shearlet_spect , 
s hear let _arg) 



for the shearlet coefficients ST. As the second argument the shearlet spectra Psi should be 
provided for faster computations, if not given, the spectra are computed with default values 
or given spectrum shearlet_spect. 

3.7 Download & Installation 

The Matlab- Version of the toolbox is available for free download at 

http : / /www . mathematik . uni-kl . de/~haeuser/FFST 

The zip-file contains all relevant files and folders. Simply unzip the archive and add the folder 
(with subfolders!) to your MATLAB-path. 

The folder FFST contains the main files for the both directions of the transform. The included 
shearlets are stored in the folder shearlets. The folder helper contains some helper functions. 
To create simple geometric structures some functions are provided in create. See contents .m 
and the comments in each file for more information. 

The following listing shows the subdirectories and the respective files 
FFST/ 

create/ 

myBall .m 
myPicture .m 
myPicture2 .m 
myRhombus . m 
mySquare .m 

helper/ 

„ conelndicator .m 

„ scalesShearsAndSpectra . m 

1 shearletScaleShear . m 

shearlets/ 

L 

bump.m 
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„ meyeraux .m 

„ meyerNewShearletSpect .m 

„ meyerScaling . m 

„ meyerScalingSpect .m 

„ meyerShearletSpect . m 

1 meyerWavelet . m 

inverseShearletTransf ormSpect .m 

shearletTransf ormSpect .m 

If everything is installed correctly run simple_example for testing. The result should look like 
Fig. 10. 



Figure 10: Result of script simple_example. 

3.8 Performance 

To evaluate the performance and the exactness of our implementation we present the following 
figures. In Fig. 11(a) we investigate the numerical tightness of the frame. The figure shows 
the difference between the square sum of the shearlets and 1, i.e., 

j'o-i 2i jo-l 

E E Efe/ + E Ei«,oi 2 + «-i 

ne{h,v} 3=0 k=-23 j=0 k=±2i 

The biggest deviation is about 8 • 10~ 15 which is 40 times the machine precision. The second 
figure Fig. 11(b) shows the difference between the original image and the back transformed 
image of a random image, i.e., the exactness of the forward and backwards transform. Here 
the biggest difference is about 2 • 10 -15 or approximately 10 times the machine precision. 
Surprisingly this is even better than the tightness of the used frame. 

Next we want to compare the speed of our implementation for different image sizes N x N 
for N = 2* with i = 5, . . . , 10. In the first run all the spectra are computed, thus, this run 
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(a) frame tightness (b) transform exactness 

Figure 11: Comparison of frame tightness and the exactness of the transform 

takes significantly longer than the second and all following runs. The different run times are 
shown in Fig. 12 with logarithmic scale. To get reasonable results we take the average time 
over 5 "first" runs (with deleting the cache, solid line) and 5 "second runs" (dashed line). The 
dotted line shows the time needed for the inverse shearlet transform. Since no spectra have to 
be computed here the time is approximately the same as for the second runs of the transform 
self. We compare the runtimes with ShearLab 2 , the only so far publicly available shearlet 
implementation. The dash-dotted line in Fig. 12 shows that ShearLab is slightly slower then 
our implementation (in the first run). Observe that no time could be measured for TV = 1024 
with ShearLab. In our implementation it is possible to compute the shearlet transform for 
arbitrary image sizes, whereas in ShearLab the transform is only possible for given image sizes. 
All tests were performed on an Intel i7 870 (Quad Core, each 2.93 GHz) with 8 GB RAM on 
Ubuntu 10.04 with Matlab R2011b (64-bit). 

3.9 Remarks 

1. In [16] and the respective implementation ShearLab a pseudo-polar Fourier transform 
is used to implement a discrete (or digital) shearlet transform. For the scale a and 
the shear s the same discretization is used. But for the translation t the authors set 
tj,k,m '■= A aj S Sjk m where we in contrast simply set t m := m (see (17)). Thus, their 
discrete shearlet becomes 

7 ( \ ?f A cT \ -2ni(u!,A a .S s . ,m) ? f a a? \ ^A a .w,m) 

i>j,k,m{u) = ip{A aj S Sjk u;)e * >* = ^{A aj S s . k u)e J 

Since the operation SJ. k A aj oj would destroy the pseudo-polar grid a "slight" adjustment 
is made and the exponential term is replaced by 

e -27ri<(0oS s - T fc )Sj jk A aj u>, m ) 

2 http://www. shearlab.org 
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Figure 12: run times for different image sizes, "first run" (solid), "second" run (dashed), inverse 
(dotted), ShearLab (dashed-dotted, no time for N = 1024). 

with 9: R \ {0} xl->KxR and 6(x, y) = (x, f ) such that 



With this adjustment the last step of the shearlet transform can be obtained with a 
standard inverse fast Fourier transform (similar as in our implementation). Unfortunately 
this is no longer related to translations of the shearlets in the time domain. 

2. We are aware of our larger oversampling factor in comparison with, e.g., ShearLab. Hav- 
ing 4 scales we obtain 61 images of the same size as the original image. But since 
shearlets are designed to detect edges in images we like to avoid any down-sampling and 
keep translation invariance. A possibility to reduce the memory usage would be to use 
the compact support of the shearlets in the frequency domain and only compute them 
on a "relevant" region. But we then would also have to store the position and size of 
each region what decreases the memory savings and would make the implementation a 
lot more complicated. 
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